This report runs base, AR1, and AR2 models (variable ~ time) for monthly SEV meteorological data using the nlme R package. Mean monthly air temperature, minimum monthly air temperature, maximum monthly air temperature, and total monthly precipitation are modeled. More sophisticated models will likely be required after this preliminary analysis. The main goal is to determine whether certain variables are changing over time for different months of the year for the SEV meteorological stations.
library(tidyverse)
library(lubridate)
library(nlme)
library(MuMIn)
library(gridExtra)
library(kableExtra)
# Load yearly data ---------------------------------------------
path_to_files <- "../data/processed_data/"
# load and only keep vars of interest
m <- read_csv(paste0(path_to_files, "met_monthly_gap_filled.csv")) %>%
mutate(sta = as.factor(sta)) %>%
select(sta:ppt)
str(m)
## tibble [1,320 × 8] (S3: tbl_df/tbl/data.frame)
## $ sta : Factor w/ 4 levels "40","42","49",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ month : num [1:1320] 1 1 1 1 1 1 1 1 1 1 ...
## $ year : num [1:1320] 1988 1989 1990 1991 1992 ...
## $ date : Date[1:1320], format: "1988-01-01" "1989-01-01" ...
## $ airt : num [1:1320] 2.291 2.291 2.291 0.615 -1.32 ...
## $ maxair: num [1:1320] 11.9 11.9 11.9 12.6 15.9 ...
## $ minair: num [1:1320] -7.2 -7.2 -7.2 -13.8 -15.1 ...
## $ ppt : num [1:1320] 0 0 3.17 7.2 16.43 ...
summary(m)
## sta month year date airt
## 40:408 Min. : 1.00 Min. :1988 Min. :1988-01-01 Min. :-1.320
## 42:396 1st Qu.: 3.75 1st Qu.:2001 1st Qu.:2001-02-22 1st Qu.: 6.362
## 49:276 Median : 6.50 Median :2008 Median :2008-03-16 Median :13.724
## 50:240 Mean : 6.50 Mean :2007 Mean :2007-07-13 Mean :13.889
## 3rd Qu.: 9.25 3rd Qu.:2015 3rd Qu.:2015-02-01 3rd Qu.:21.895
## Max. :12.00 Max. :2021 Max. :2021-12-01 Max. :29.569
## maxair minair ppt
## Min. : 8.985 Min. :-29.930 Min. : 0.000
## 1st Qu.:21.475 1st Qu.: -9.523 1st Qu.: 4.104
## Median :28.725 Median : -2.680 Median : 13.524
## Mean :28.059 Mean : -1.336 Mean : 21.875
## 3rd Qu.:34.655 3rd Qu.: 7.400 3rd Qu.: 31.050
## Max. :48.720 Max. : 17.910 Max. :183.147
Split data into individual met stations and add a time variable for nlme modeling.
# function to select station and variable and add 'time' var for a specific month of the year
prepare_data_for_nlme <- function(data, sta, month) {
# data = dataframe to use
# sta = met station (in quotes)
# month = month of year (in quotes)
data %>%
filter(sta == {{ sta }} & month == {{ month }}) %>%
arrange(date) %>%
mutate(time = as.numeric(as.factor(year)))
}
This code produces several functions that will be useful for producing cleaner code.
# function for plotting preliminary monthly graph -----------------------------------------------------------
plot_monthly_prelim <- function(data, var, title, y_axis_label) {
# data = dataset
# var = variable (not in quotes)
# title = title for graph (in quotes)
# y_axis_label = label for y-axis (in quotes)
ggplot(data, aes(x = date, y = {{ var }}, col = sta)) +
geom_line(color = "burlywood") +
geom_smooth(method = "lm", color = "black", size = 0.6) +
facet_wrap(~ sta) +
labs(title = title,
x = "Month",
y = y_axis_label)
}
by_month_lm <- function(data, var, station, title, y_axis_label) {
# data = dataset
# var = variable to plot (not in quotes)
# station = met station (in quotes)
# title = graph title
# y_axis_label = label for the graph's y-axis
m %>%
filter(sta == station) %>%
ggplot(., aes(x = year, y = {{ var }})) +
geom_line() +
geom_smooth(method = "lm") +
facet_wrap(~ month, scales = "free_y") +
ggtitle({{ title }})
}
# function for plotting graphs with best nlme model results ------------------------------------------------
plot_monthly_results <- function(data, var, coeffs, title, y_axis_label) {
# data = dataset
# var = variable (not in quotes)
# coeffs = object containing coefficients and p-values of best model
# title = title for graph (in quotes)
# y_axis_label = label for y-axis (in quotes)
ggplot(data, aes(x = date, y = {{ var }})) +
geom_line(color = "burlywood") +
geom_smooth(method = "lm", color = "black", size = 0.6) +
labs(title = title,
x = "Month",
y = y_axis_label,
caption = paste("slope = ", round(coeffs[2,1], 3), "\n(p-value = ", coeffs[2,2], ")"))
}
# functions to run nlme models ----------------------------------------------------------
# base model function
base_nlme_model <- function(data, var, month) {
# data = dataframe
# var = var to model
# month = month of year
data <- data %>% filter(month == month)
model_specification = as.formula(paste0(var, " ~ time"))
gls(model = model_specification,
data = data,
method = "ML",
na.action = na.omit)
}
# AR1 model function
AR1_nlme_model <- function(data, var, month) {
# data = dataframe
# var = var to model
# month = month of year
data <- data %>% filter(month == month)
model_specification = as.formula(paste0(var, " ~ time"))
gls(model = model_specification,
data = data,
correlation = corARMA(form = ~ time, p = 1),
method = "ML",
na.action = na.omit)
}
# AR2 model function
AR2_nlme_model <- function(data, var, month) {
# data = dataframe
# var = var to model
# month = month of year
data <- data %>% filter(month == month)
model_specification = as.formula(paste0(var, " ~ time"))
gls(model = model_specification,
data = data,
correlation = corARMA(form = ~ time, p = 2),
method = "ML",
na.action = na.omit)
}
plot_monthly_prelim(m, airt, "Monthly Mean Air Temperature", "Temperature (C)")
## `geom_smooth()` using formula 'y ~ x'
by_month_lm(m, airt, "40", "Met 40 - Deep Well \nMean Monthly Air Temperature", "Temperature (C)")
## `geom_smooth()` using formula 'y ~ x'
m40_6 <- prepare_data_for_nlme(m, "40", "6")
m40_6_mbase <- base_nlme_model(m40_6, "airt", 6)
m40_6_mAR1 <- AR1_nlme_model(m40_6, "airt", 6)
m40_5_mAR2 <- AR2_nlme_model(m40_6, "airt", 6)
AICc scores:
AICc(m40_6_mbase, m40_6_mAR1, m40_5_mAR2) %>%
kbl() %>%
kable_styling(full_width = FALSE, position = "left")
| df | AICc | |
|---|---|---|
| m40_6_mbase | 3 | 110.5176 |
| m40_6_mAR1 | 4 | 112.6553 |
| m40_5_mAR2 | 5 | 114.8095 |
summary(m40_6_mbase)
## Generalized least squares fit by maximum likelihood
## Model: model_specification
## Data: data
## AIC BIC logLik
## 109.7176 114.2966 -51.85878
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 22.828769 0.4020514 56.78072 0.0000
## time 0.059215 0.0200400 2.95484 0.0058
##
## Correlation:
## (Intr)
## time -0.872
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -1.8703743 -0.7313602 0.2242368 0.6785338 1.9307997
##
## Residual standard error: 1.112177
## Degrees of freedom: 34 total; 32 residual
m40_6_coeffs <- summary(m40_6_mbase)$tTable[,c(1, 4)]
(m40_6_plot <- plot_monthly_results(m40_6, airt, m40_6_coeffs, "Met 40 - Deep Well \nJune Mean Air Temperature", "Temperature (C)"))
## `geom_smooth()` using formula 'y ~ x'
by_month_lm(m, airt, "42", "Met 42 - Deep Well \nMean Monthly Air Temperature", "Temperature (C)")
## `geom_smooth()` using formula 'y ~ x'
m42_6 <- prepare_data_for_nlme(m, "42", "6")
m42_6_mbase <- base_nlme_model(m42_6, "airt", 6)
m42_6_mAR1 <- AR1_nlme_model(m42_6, "airt", 6)
m42_5_mAR2 <- AR2_nlme_model(m42_6, "airt", 6)
AICc scores:
AICc(m42_6_mbase, m42_6_mAR1, m42_5_mAR2) %>%
kbl() %>%
kable_styling(full_width = FALSE, position = "left")
| df | AICc | |
|---|---|---|
| m42_6_mbase | 3 | 105.9467 |
| m42_6_mAR1 | 4 | 108.3905 |
| m42_5_mAR2 | 5 | 111.1827 |
summary(m42_6_mbase)
## Generalized least squares fit by maximum likelihood
## Model: model_specification
## Data: data
## AIC BIC logLik
## 105.1191 109.6086 -49.55955
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 21.896715 0.3992867 54.83958 0.0000
## time 0.046355 0.0204920 2.26211 0.0308
##
## Correlation:
## (Intr)
## time -0.872
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -2.10520508 -0.68571065 -0.05950743 0.72953403 2.59403949
##
## Residual standard error: 1.086396
## Degrees of freedom: 33 total; 31 residual
m42_6_coeffs <- summary(m42_6_mbase)$tTable[,c(1, 4)]
(m42_6_plot <- plot_monthly_results(m42_6, airt, m42_6_coeffs, "Met 42 - Cerro Montoso \nJune Mean Air Temperature", "Temperature (C)"))
## `geom_smooth()` using formula 'y ~ x'
by_month_lm(m49_6, airt, "49", "Met 49 - Deep Well \nMean Monthly Air Temperature", "Temperature (C)")
## `geom_smooth()` using formula 'y ~ x'
m49_6 <- prepare_data_for_nlme(m, "49", "6")
m49_6_mbase <- base_nlme_model(m49_6, "airt", 6)
m49_6_mAR1 <- AR1_nlme_model(m49_6, "airt", 6)
m49_5_mAR2 <- AR2_nlme_model(m49_6, "airt", 6)
AICc scores:
AICc(m49_6_mbase, m49_6_mAR1, m49_5_mAR2) %>%
kbl() %>%
kable_styling(full_width = FALSE, position = "left")
| df | AICc | |
|---|---|---|
| m49_6_mbase | 3 | 71.53132 |
| m49_6_mAR1 | 4 | 74.26520 |
| m49_5_mAR2 | 5 | 77.55554 |
summary(m49_6_mbase)
## Generalized least squares fit by maximum likelihood
## Model: model_specification
## Data: data
## AIC BIC logLik
## 70.26816 73.67465 -32.13408
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 23.986467 0.4413401 54.34917 0.0000
## time 0.066771 0.0321880 2.07440 0.0505
##
## Correlation:
## (Intr)
## time -0.875
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -2.2133523 -0.7132157 0.1294271 0.5771409 1.7943170
##
## Residual standard error: 0.9784315
## Degrees of freedom: 23 total; 21 residual
m49_6_coeffs <- summary(m49_6_mbase)$tTable[,c(1, 4)]
(m49_6_plot <- plot_monthly_results(m49_6, airt, m49_6_coeffs, "Met 49 - Five Points \nJune Mean Air Temperature", "Temperature (C)"))
## `geom_smooth()` using formula 'y ~ x'
by_month_lm(m, airt, "50", "Met 50 - Deep Well \nMean Monthly Air Temperature", "Temperature (C)")
## `geom_smooth()` using formula 'y ~ x'
m50_6 <- prepare_data_for_nlme(m, "50", "6")
m50_6_mbase <- base_nlme_model(m50_6, "airt", 6)
m50_6_mAR1 <- AR1_nlme_model(m50_6, "airt", 6)
m50_5_mAR2 <- AR2_nlme_model(m50_6, "airt", 6)
AICc scores:
AICc(m50_6_mbase, m50_6_mAR1, m50_5_mAR2) %>%
kbl() %>%
kable_styling(full_width = FALSE, position = "left")
| df | AICc | |
|---|---|---|
| m50_6_mbase | 3 | 68.93686 |
| m50_6_mAR1 | 4 | 71.45522 |
| m50_5_mAR2 | 5 | 75.03545 |
summary(m50_6_mbase)
## Generalized least squares fit by maximum likelihood
## Model: model_specification
## Data: data
## AIC BIC logLik
## 67.43686 70.42406 -30.71843
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 25.066671 0.5504260 45.54049 0.0000
## time 0.072846 0.0459487 1.58538 0.1303
##
## Correlation:
## (Intr)
## time -0.877
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -2.30015742 -0.66646539 0.08760363 0.68926206 1.49495780
##
## Residual standard error: 1.1241
## Degrees of freedom: 20 total; 18 residual
m50_6_coeffs <- summary(m50_6_mbase)$tTable[,c(1, 4)]
(m50_6_plot <- plot_monthly_results(m50_6, airt, m50_6_coeffs, "Met 50 - Blue Grama \nJune Mean Air Temperature", "Temperature (C)"))
## `geom_smooth()` using formula 'y ~ x'
# airt results
(airt_6_graphs <- grid.arrange(m40_6_plot, m42_6_plot, m49_6_plot, m50_6_plot, ncol=2))
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
## TableGrob (2 x 2) "arrange": 4 grobs
## z cells name grob
## 1 1 (1-1,1-1) arrange gtable[layout]
## 2 2 (1-1,2-2) arrange gtable[layout]
## 3 3 (2-2,1-1) arrange gtable[layout]
## 4 4 (2-2,2-2) arrange gtable[layout]
plot_monthly_prelim(m, maxair, "Monthly Maximum Air Temperature", "Temperature (C)")
## `geom_smooth()` using formula 'y ~ x'
by_month_lm(m, maxair, "40", "Met 40 - Deep Well \nMaximum Monthly Air Temperature", "Temperature (C)")
## `geom_smooth()` using formula 'y ~ x'
m40_6 <- prepare_data_for_nlme(m, "40", "6")
m40_6_mbase <- base_nlme_model(m40_6, "maxair", 6)
m40_6_mAR1 <- AR1_nlme_model(m40_6, "maxair", 6)
m40_5_mAR2 <- AR2_nlme_model(m40_6, "maxair", 6)
AICc scores:
AICc(m40_6_mbase, m40_6_mAR1, m40_5_mAR2) %>%
kbl() %>%
kable_styling(full_width = FALSE, position = "left")
| df | AICc | |
|---|---|---|
| m40_6_mbase | 3 | 131.0670 |
| m40_6_mAR1 | 4 | 132.3726 |
| m40_5_mAR2 | 5 | 134.7689 |
summary(m40_6_mbase)
## Generalized least squares fit by maximum likelihood
## Model: model_specification
## Data: data
## AIC BIC logLik
## 130.267 134.8461 -62.13351
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 36.74887 0.5439069 67.56463 0.0000
## time 0.07725 0.0271108 2.84947 0.0076
##
## Correlation:
## (Intr)
## time -0.872
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -2.35174301 -0.43438926 0.07110273 0.48204954 2.65213862
##
## Residual standard error: 1.504586
## Degrees of freedom: 34 total; 32 residual
m40_6_coeffs <- summary(m40_6_mbase)$tTable[,c(1, 4)]
(m40_6_plot <- plot_monthly_results(m40_6, maxair, m40_6_coeffs, "Met 40 - Deep Well \nJune Maximum Air Temperature", "Temperature (C)"))
## `geom_smooth()` using formula 'y ~ x'
by_month_lm(m, maxair, "42", "Met 42 - Deep Well \nMaximum Monthly Air Temperature", "Temperature (C)")
## `geom_smooth()` using formula 'y ~ x'
m42_6 <- prepare_data_for_nlme(m, "42", "6")
m42_6_mbase <- base_nlme_model(m42_6, "maxair", 6)
m42_6_mAR1 <- AR1_nlme_model(m42_6, "maxair", 6)
m42_5_mAR2 <- AR2_nlme_model(m42_6, "maxair", 6)
AICc scores:
AICc(m42_6_mbase, m42_6_mAR1, m42_5_mAR2) %>%
kbl() %>%
kable_styling(full_width = FALSE, position = "left")
| df | AICc | |
|---|---|---|
| m42_6_mbase | 3 | 130.7457 |
| m42_6_mAR1 | 4 | 133.1964 |
| m42_5_mAR2 | 5 | 135.4923 |
summary(m42_6_mbase)
## Generalized least squares fit by maximum likelihood
## Model: model_specification
## Data: data
## AIC BIC logLik
## 129.9181 134.4076 -61.95906
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 34.12709 0.5813903 58.69911 0.0000
## time 0.08010 0.0298378 2.68441 0.0116
##
## Correlation:
## (Intr)
## time -0.872
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -1.76318154 -0.69378514 -0.02516334 0.53254899 3.12435430
##
## Residual standard error: 1.581872
## Degrees of freedom: 33 total; 31 residual
m42_6_coeffs <- summary(m42_6_mbase)$tTable[,c(1, 4)]
(m42_6_plot <- plot_monthly_results(m42_6, maxair, m42_6_coeffs, "Met 42 - Cerro Montoso \nJune Maximum Air Temperature", "Temperature (C)"))
## `geom_smooth()` using formula 'y ~ x'
by_month_lm(m49_6, maxair, "49", "Met 49 - Deep Well \nMaximum Monthly Air Temperature", "Temperature (C)")
## `geom_smooth()` using formula 'y ~ x'
m49_6 <- prepare_data_for_nlme(m, "49", "6")
m49_6_mbase <- base_nlme_model(m49_6, "maxair", 6)
m49_6_mAR1 <- AR1_nlme_model(m49_6, "maxair", 6)
m49_5_mAR2 <- AR2_nlme_model(m49_6, "maxair", 6)
AICc scores:
AICc(m49_6_mbase, m49_6_mAR1, m49_5_mAR2) %>%
kbl() %>%
kable_styling(full_width = FALSE, position = "left")
| df | AICc | |
|---|---|---|
| m49_6_mbase | 3 | 82.12418 |
| m49_6_mAR1 | 4 | 84.65632 |
| m49_5_mAR2 | 5 | 87.96102 |
summary(m49_6_mbase)
## Generalized least squares fit by maximum likelihood
## Model: model_specification
## Data: data
## AIC BIC logLik
## 80.86102 84.26751 -37.43051
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 37.52901 0.5556259 67.54367 0.0000
## time 0.05914 0.0405232 1.45942 0.1592
##
## Correlation:
## (Intr)
## time -0.875
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -2.5406398 -0.3613768 0.2014793 0.6971315 1.8703417
##
## Residual standard error: 1.231798
## Degrees of freedom: 23 total; 21 residual
m49_6_coeffs <- summary(m49_6_mbase)$tTable[,c(1, 4)]
(m49_6_plot <- plot_monthly_results(m49_6, maxair, m49_6_coeffs, "Met 49 - Five Points \nJune Maximum Air Temperature", "Temperature (C)"))
## `geom_smooth()` using formula 'y ~ x'
by_month_lm(m, maxair, "50", "Met 50 - Deep Well \nMaximum Monthly Air Temperature", "Temperature (C)")
## `geom_smooth()` using formula 'y ~ x'
m50_6 <- prepare_data_for_nlme(m, "50", "6")
m50_6_mbase <- base_nlme_model(m50_6, "maxair", 6)
m50_6_mAR1 <- AR1_nlme_model(m50_6, "maxair", 6)
m50_5_mAR2 <- AR2_nlme_model(m50_6, "maxair", 6)
AICc scores:
AICc(m50_6_mbase, m50_6_mAR1, m50_5_mAR2) %>%
kbl() %>%
kable_styling(full_width = FALSE, position = "left")
| df | AICc | |
|---|---|---|
| m50_6_mbase | 3 | 78.64979 |
| m50_6_mAR1 | 4 | 81.79530 |
| m50_5_mAR2 | 5 | 85.41433 |
summary(m50_6_mbase)
## Generalized least squares fit by maximum likelihood
## Model: model_specification
## Data: data
## AIC BIC logLik
## 77.14979 80.13699 -35.57489
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 37.54153 0.7017069 53.50030 0.0000
## time 0.12119 0.0585774 2.06885 0.0532
##
## Correlation:
## (Intr)
## time -0.877
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -2.3244313 -0.5513157 0.1981575 0.5670821 1.9498376
##
## Residual standard error: 1.433052
## Degrees of freedom: 20 total; 18 residual
m50_6_coeffs <- summary(m50_6_mbase)$tTable[,c(1, 4)]
(m50_6_plot <- plot_monthly_results(m50_6, maxair, m50_6_coeffs, "Met 50 - Blue Grama \nJune Maximum Air Temperature", "Temperature (C)"))
## `geom_smooth()` using formula 'y ~ x'
# maxair results
(maxair_6_graphs <- grid.arrange(m40_6_plot, m42_6_plot, m49_6_plot, m50_6_plot, ncol=2))
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
## TableGrob (2 x 2) "arrange": 4 grobs
## z cells name grob
## 1 1 (1-1,1-1) arrange gtable[layout]
## 2 2 (1-1,2-2) arrange gtable[layout]
## 3 3 (2-2,1-1) arrange gtable[layout]
## 4 4 (2-2,2-2) arrange gtable[layout]
plot_monthly_prelim(m, minair, "Monthly Minimum Air Temperature", "Temperature (C)")
## `geom_smooth()` using formula 'y ~ x'
by_month_lm(m, minair, "40", "Met 40 - Deep Well \nMinimum Monthly Air Temperature", "Temperature (C)")
## `geom_smooth()` using formula 'y ~ x'
m40_6 <- prepare_data_for_nlme(m, "40", "6")
m40_6_mbase <- base_nlme_model(m40_6, "minair", 6)
m40_6_mAR1 <- AR1_nlme_model(m40_6, "minair", 6)
m40_5_mAR2 <- AR2_nlme_model(m40_6, "minair", 6)
AICc scores:
AICc(m40_6_mbase, m40_6_mAR1, m40_5_mAR2) %>%
kbl() %>%
kable_styling(full_width = FALSE, position = "left")
| df | AICc | |
|---|---|---|
| m40_6_mbase | 3 | 168.6741 |
| m40_6_mAR1 | 4 | 171.2504 |
| m40_5_mAR2 | 5 | 173.8996 |
summary(m40_6_mbase)
## Generalized least squares fit by maximum likelihood
## Model: model_specification
## Data: data
## AIC BIC logLik
## 167.8741 172.4532 -80.93704
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 7.833538 0.9456032 8.284170 0.0000
## time -0.029816 0.0471331 -0.632596 0.5315
##
## Correlation:
## (Intr)
## time -0.872
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -2.25544260 -0.73370381 -0.09415726 0.75553571 1.62187112
##
## Residual standard error: 2.615781
## Degrees of freedom: 34 total; 32 residual
m40_6_coeffs <- summary(m40_6_mbase)$tTable[,c(1, 4)]
(m40_6_plot <- plot_monthly_results(m40_6, minair, m40_6_coeffs, "Met 40 - Deep Well \nJune Minimum Air Temperature", "Temperature (C)"))
## `geom_smooth()` using formula 'y ~ x'
by_month_lm(m, minair, "42", "Met 42 - Deep Well \nMinimum Monthly Air Temperature", "Temperature (C)")
## `geom_smooth()` using formula 'y ~ x'
m42_6 <- prepare_data_for_nlme(m, "42", "6")
m42_6_mbase <- base_nlme_model(m42_6, "minair", 6)
m42_6_mAR1 <- AR1_nlme_model(m42_6, "minair", 6)
m42_5_mAR2 <- AR2_nlme_model(m42_6, "minair", 6)
AICc scores:
AICc(m42_6_mbase, m42_6_mAR1, m42_5_mAR2) %>%
kbl() %>%
kable_styling(full_width = FALSE, position = "left")
| df | AICc | |
|---|---|---|
| m42_6_mbase | 3 | 153.6229 |
| m42_6_mAR1 | 4 | 155.3432 |
| m42_5_mAR2 | 5 | 155.8929 |
summary(m42_6_mbase)
## Generalized least squares fit by maximum likelihood
## Model: model_specification
## Data: data
## AIC BIC logLik
## 152.7953 157.2848 -73.39764
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 8.669623 0.8222512 10.54376 0.0000
## time -0.039671 0.0421991 -0.94009 0.3544
##
## Correlation:
## (Intr)
## time -0.872
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -2.2567953 -0.6711610 -0.1724796 0.7295458 2.1228551
##
## Residual standard error: 2.237216
## Degrees of freedom: 33 total; 31 residual
m42_6_coeffs <- summary(m42_6_mbase)$tTable[,c(1, 4)]
(m42_6_plot <- plot_monthly_results(m42_6, minair, m42_6_coeffs, "Met 42 - Cerro Montoso \nJune Minimum Air Temperature", "Temperature (C)"))
## `geom_smooth()` using formula 'y ~ x'
by_month_lm(m49_6, minair, "49", "Met 49 - Deep Well \nMinimum Monthly Air Temperature", "Temperature (C)")
## `geom_smooth()` using formula 'y ~ x'
m49_6 <- prepare_data_for_nlme(m, "49", "6")
m49_6_mbase <- base_nlme_model(m49_6, "minair", 6)
m49_6_mAR1 <- AR1_nlme_model(m49_6, "minair", 6)
m49_5_mAR2 <- AR2_nlme_model(m49_6, "minair", 6)
AICc scores:
AICc(m49_6_mbase, m49_6_mAR1, m49_5_mAR2) %>%
kbl() %>%
kable_styling(full_width = FALSE, position = "left")
| df | AICc | |
|---|---|---|
| m49_6_mbase | 3 | 103.5652 |
| m49_6_mAR1 | 4 | 101.6278 |
| m49_5_mAR2 | 5 | 104.2407 |
summary(m49_6_mAR1)
## Generalized least squares fit by maximum likelihood
## Model: model_specification
## Data: data
## AIC BIC logLik
## 99.40558 103.9476 -45.70279
##
## Correlation Structure: AR(1)
## Formula: ~time
## Parameter estimate(s):
## Phi
## -0.4759774
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 8.610738 0.5561496 15.482773 0.0000
## time 0.058075 0.0408268 1.422464 0.1696
##
## Correlation:
## (Intr)
## time -0.881
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -1.97613576 -0.57321163 -0.02887857 0.77814343 1.67707143
##
## Residual standard error: 1.995719
## Degrees of freedom: 23 total; 21 residual
m49_6_coeffs <- summary(m49_6_mAR1)$tTable[,c(1, 4)]
(m49_6_plot <- plot_monthly_results(m49_6, minair, m49_6_coeffs, "Met 49 - Five Points \nJune Minimum Air Temperature", "Temperature (C)"))
## `geom_smooth()` using formula 'y ~ x'
by_month_lm(m, minair, "50", "Met 50 - Deep Well \nMinimum Monthly Air Temperature", "Temperature (C)")
## `geom_smooth()` using formula 'y ~ x'
m50_6 <- prepare_data_for_nlme(m, "50", "6")
m50_6_mbase <- base_nlme_model(m50_6, "minair", 6)
m50_6_mAR1 <- AR1_nlme_model(m50_6, "minair", 6)
m50_5_mAR2 <- AR2_nlme_model(m50_6, "minair", 6)
AICc scores:
AICc(m50_6_mbase, m50_6_mAR1, m50_5_mAR2) %>%
kbl() %>%
kable_styling(full_width = FALSE, position = "left")
| df | AICc | |
|---|---|---|
| m50_6_mbase | 3 | 90.80023 |
| m50_6_mAR1 | 4 | 90.42714 |
| m50_5_mAR2 | 5 | 93.73701 |
summary(m50_6_mAR1)
## Generalized least squares fit by maximum likelihood
## Model: model_specification
## Data: data
## AIC BIC logLik
## 87.76048 91.74341 -39.88024
##
## Correlation Structure: AR(1)
## Formula: ~time
## Parameter estimate(s):
## Phi
## 0.4078502
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 12.354002 1.3459652 9.178545 0.000
## time -0.104734 0.1105605 -0.947298 0.356
##
## Correlation:
## (Intr)
## time -0.862
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -1.9317235 -0.5736404 0.1539656 0.5339572 2.3574953
##
## Residual standard error: 1.937679
## Degrees of freedom: 20 total; 18 residual
m50_6_coeffs <- summary(m50_6_mAR1)$tTable[,c(1, 4)]
(m50_6_plot <- plot_monthly_results(m50_6, minair, m50_6_coeffs, "Met 50 - Blue Grama \nJune Minimum Air Temperature", "Temperature (C)"))
## `geom_smooth()` using formula 'y ~ x'
# minair results
(minair_6_graphs <- grid.arrange(m40_6_plot, m42_6_plot, m49_6_plot, m50_6_plot, ncol=2))
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
## TableGrob (2 x 2) "arrange": 4 grobs
## z cells name grob
## 1 1 (1-1,1-1) arrange gtable[layout]
## 2 2 (1-1,2-2) arrange gtable[layout]
## 3 3 (2-2,1-1) arrange gtable[layout]
## 4 4 (2-2,2-2) arrange gtable[layout]
minair_6_graphs
## TableGrob (2 x 2) "arrange": 4 grobs
## z cells name grob
## 1 1 (1-1,1-1) arrange gtable[layout]
## 2 2 (1-1,2-2) arrange gtable[layout]
## 3 3 (2-2,1-1) arrange gtable[layout]
## 4 4 (2-2,2-2) arrange gtable[layout]
airt_6_graphs
## TableGrob (2 x 2) "arrange": 4 grobs
## z cells name grob
## 1 1 (1-1,1-1) arrange gtable[layout]
## 2 2 (1-1,2-2) arrange gtable[layout]
## 3 3 (2-2,1-1) arrange gtable[layout]
## 4 4 (2-2,2-2) arrange gtable[layout]
maxair_6_graphs
## TableGrob (2 x 2) "arrange": 4 grobs
## z cells name grob
## 1 1 (1-1,1-1) arrange gtable[layout]
## 2 2 (1-1,2-2) arrange gtable[layout]
## 3 3 (2-2,1-1) arrange gtable[layout]
## 4 4 (2-2,2-2) arrange gtable[layout]
plot_monthly_prelim(m, ppt, "Monthly Precipitation", "Precipitation (mm)")
## `geom_smooth()` using formula 'y ~ x'
by_month_lm(m, ppt, "40", "Met 40 - Deep Well \n Monthly Precipitation", "Precipitation (mm)")
## `geom_smooth()` using formula 'y ~ x'
m40_6 <- prepare_data_for_nlme(m, "40", "6")
m40_6_mbase <- base_nlme_model(m40_6, "ppt", 6)
m40_6_mAR1 <- AR1_nlme_model(m40_6, "ppt", 6)
m40_5_mAR2 <- AR2_nlme_model(m40_6, "ppt", 6)
AICc scores:
AICc(m40_6_mbase, m40_6_mAR1, m40_5_mAR2) %>%
kbl() %>%
kable_styling(full_width = FALSE, position = "left")
| df | AICc | |
|---|---|---|
| m40_6_mbase | 3 | 303.5660 |
| m40_6_mAR1 | 4 | 306.0801 |
| m40_5_mAR2 | 5 | 306.6254 |
summary(m40_6_mbase)
## Generalized least squares fit by maximum likelihood
## Model: model_specification
## Data: data
## AIC BIC logLik
## 302.766 307.3451 -148.383
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 21.708889 6.874179 3.1580337 0.0035
## time -0.282118 0.342640 -0.8233652 0.4164
##
## Correlation:
## (Intr)
## time -0.872
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -1.1119553 -0.7733032 -0.3014592 0.3845706 2.9864295
##
## Residual standard error: 19.01574
## Degrees of freedom: 34 total; 32 residual
m40_6_coeffs <- summary(m40_6_mbase)$tTable[,c(1, 4)]
(m40_6_plot <- plot_monthly_results(m40_6, ppt, m40_6_coeffs, "Met 40 - Deep Well \nJune Precipitation", "Precipitation (mm)"))
## `geom_smooth()` using formula 'y ~ x'
by_month_lm(m, ppt, "42", "Met 42 - Deep Well \n Monthly Precipitation", "Precipitation (mm)")
## `geom_smooth()` using formula 'y ~ x'
m42_6 <- prepare_data_for_nlme(m, "42", "6")
m42_6_mbase <- base_nlme_model(m42_6, "ppt", 6)
m42_6_mAR1 <- AR1_nlme_model(m42_6, "ppt", 6)
m42_5_mAR2 <- AR2_nlme_model(m42_6, "ppt", 6)
AICc scores:
AICc(m42_6_mbase, m42_6_mAR1, m42_5_mAR2) %>%
kbl() %>%
kable_styling(full_width = FALSE, position = "left")
| df | AICc | |
|---|---|---|
| m42_6_mbase | 3 | 303.7262 |
| m42_6_mAR1 | 4 | 306.0666 |
| m42_5_mAR2 | 5 | 305.4861 |
summary(m42_6_mbase)
## Generalized least squares fit by maximum likelihood
## Model: model_specification
## Data: data
## AIC BIC logLik
## 302.8986 307.3881 -148.4493
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 26.517472 7.993137 3.317530 0.0023
## time -0.355085 0.410219 -0.865599 0.3934
##
## Correlation:
## (Intr)
## time -0.872
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -1.2029751 -0.6962164 -0.3386889 0.3329534 3.9586599
##
## Residual standard error: 21.74807
## Degrees of freedom: 33 total; 31 residual
m42_6_coeffs <- summary(m42_6_mbase)$tTable[,c(1, 4)]
(m42_6_plot <- plot_monthly_results(m42_6, ppt, m42_6_coeffs, "Met 42 - Cerro Montoso \nJune Precipitation", "Precipitation (mm)"))
## `geom_smooth()` using formula 'y ~ x'
by_month_lm(m49_6, ppt, "49", "Met 49 - Deep Well \n Monthly Precipitation", "Precipitation (mm)")
## `geom_smooth()` using formula 'y ~ x'
m49_6 <- prepare_data_for_nlme(m, "49", "6")
m49_6_mbase <- base_nlme_model(m49_6, "ppt", 6)
m49_6_mAR1 <- AR1_nlme_model(m49_6, "ppt", 6)
m49_5_mAR2 <- AR2_nlme_model(m49_6, "ppt", 6)
AICc scores:
AICc(m49_6_mbase, m49_6_mAR1, m49_5_mAR2) %>%
kbl() %>%
kable_styling(full_width = FALSE, position = "left")
| df | AICc | |
|---|---|---|
| m49_6_mbase | 3 | 177.3561 |
| m49_6_mAR1 | 4 | 178.3878 |
| m49_5_mAR2 | 5 | 175.5247 |
summary(m49_6_mbase)
## Generalized least squares fit by maximum likelihood
## Model: model_specification
## Data: data
## AIC BIC logLik
## 176.0929 179.4994 -85.04645
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 11.236151 4.404374 2.5511347 0.0186
## time 0.111922 0.321222 0.3484245 0.7310
##
## Correlation:
## (Intr)
## time -0.875
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -1.2997479 -0.8456391 -0.3466543 0.7219178 2.1507099
##
## Residual standard error: 9.764303
## Degrees of freedom: 23 total; 21 residual
m49_6_coeffs <- summary(m49_6_mbase)$tTable[,c(1, 4)]
(m49_6_plot <- plot_monthly_results(m49_6, ppt, m49_6_coeffs, "Met 49 - Five Points \nJune Precipitation", "Precipitation (mm)"))
## `geom_smooth()` using formula 'y ~ x'
by_month_lm(m, ppt, "50", "Met 50 - Deep Well \n Monthly Precipitation", "Precipitation (mm)")
## `geom_smooth()` using formula 'y ~ x'
m50_6 <- prepare_data_for_nlme(m, "50", "6")
m50_6_mbase <- base_nlme_model(m50_6, "ppt", 6)
m50_6_mAR1 <- AR1_nlme_model(m50_6, "ppt", 6)
m50_5_mAR2 <- AR2_nlme_model(m50_6, "ppt", 6)
AICc scores:
AICc(m50_6_mbase, m50_6_mAR1, m50_5_mAR2) %>%
kbl() %>%
kable_styling(full_width = FALSE, position = "left")
| df | AICc | |
|---|---|---|
| m50_6_mbase | 3 | 176.3212 |
| m50_6_mAR1 | 4 | 177.2593 |
| m50_5_mAR2 | 5 | 178.6292 |
summary(m50_6_mbase)
## Generalized least squares fit by maximum likelihood
## Model: model_specification
## Data: data
## AIC BIC logLik
## 174.8212 177.8084 -84.41058
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 21.677368 8.065088 2.6878030 0.0150
## time -0.611654 0.673261 -0.9084954 0.3756
##
## Correlation:
## (Intr)
## time -0.877
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -1.0704234 -0.5841853 -0.2183698 0.4232182 3.4459480
##
## Residual standard error: 16.47082
## Degrees of freedom: 20 total; 18 residual
m50_6_coeffs <- summary(m50_6_mbase)$tTable[,c(1, 4)]
(m50_6_plot <- plot_monthly_results(m50_6, ppt, m50_6_coeffs, "Met 50 - Blue Grama \nJune Precipitation", "Precipitation (mm)"))
## `geom_smooth()` using formula 'y ~ x'
# ppt results
(ppt_6_graphs <- grid.arrange(m40_6_plot, m42_6_plot, m49_6_plot, m50_6_plot, ncol=2))
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
## TableGrob (2 x 2) "arrange": 4 grobs
## z cells name grob
## 1 1 (1-1,1-1) arrange gtable[layout]
## 2 2 (1-1,2-2) arrange gtable[layout]
## 3 3 (2-2,1-1) arrange gtable[layout]
## 4 4 (2-2,2-2) arrange gtable[layout]